home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / The World of Computer Software.iso / dec92.zip / 1012090A < prev    next >
Text File  |  1992-10-13  |  4KB  |  124 lines
  1. /* Copyright (c) 1985 Regents of the University of California.
  2.  
  3.   Use and reproduction of this software are granted in accordance with
  4.   the terms and conditions specified in the Berkeley Software License
  5.   Agreement (in particular, this entails acknowledgement of the
  6.   programs' source, and inclusion of this notice) with the additional
  7.   understanding that all recipients should regard themselves as
  8.   participants in an ongoing research project and hence should feel
  9.   obligated to report their experiences (good or bad) with these
  10.   elementary function codes, using "sendbug 4bsd-bugs@BERKELEY", to the
  11.   authors.
  12.  
  13.  */
  14.  
  15. #ifndef lint
  16. static char sccsid[] = "@(#)atan2.c     (T Prince) 1/3/92";
  17. #endif                          /* not lint */
  18.  
  19. /* ATAN2(Y,X) RETURN ARG (X+iY)
  20.  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  21.  * CODED IN C BY K.C. NG, 1/8/85;
  22.  * REVISED BY K.C. NG on 2/7/85, 2/13/85, 3/7/85, 3/30/85, 6/29/85.
  23.  * simplified by T C Prince 12/91
  24.  *
  25.  * Method :
  26.  1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
  27.  2. Reduce x to positive by (if x and y are unexceptional):
  28.                 ARG (x+iy) = arctan(y/x)           ... if x > 0,
  29.                 ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
  30.  3.  According to the integer k=4t+0.25 truncated, t=y/x, the argument
  31.         is further reduced to one of the following intervals and the
  32.         arctangent of y/x is evaluated by the corresponding formula:
  33.  
  34.         [0,7/16]        atan(y/x) = t - t^3*(a1+t^2*(a2+...a10+t^2*a11)
  35.         [7/16,11/16]    atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) )
  36.         [11/16.19/16]   atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) )
  37.         [19/16,39/16]   atan(y/x) = atan(3/2) + atan((y-1.5x)/(x+1.5y))
  38.         [39/16,INF]     atan(y/x) = atan(INF) + atan( -x/y )
  39.  
  40.   The decimal values may be used, provided that the compiler will
  41.   convert from decimal to binary accurately enough.
  42.  
  43.  */
  44.  
  45. const static double
  46.  athfhi = .46364760900080609352, athflo = 4.6249969567426939759E-18,
  47.  at1fhi = .98279372324732905408, at1flo = -2.4407677060164810007E-17,
  48.  PI = 3.14159265358979323846,
  49. /* Chebyshev economized coefficients, 56 bit precision */
  50.  a1 = .33333333333332713, a2 = -.19999999999844163,
  51.  a3 = .14285714270355533, a4 = -.11111110327242694,
  52.  a5 = .09090885408512523, a6 = -.0769185192745614,
  53.  a7 = .06660850184641357, a8 = -.05832239923826337,
  54.  a9 = .04971908607172078, a10 = -.03642599200325829,
  55.  a11 = .01618847031840557;
  56.  
  57. #include <math.h>
  58.  
  59. double (atan) (x)
  60. double x;
  61. {
  62.     return atan2(x, 1.);
  63. }
  64.  
  65. double atan2(y, x)
  66. double y, x;
  67. {
  68.     double t, z, signy, hi, lo, z2;
  69.     int k, m, signx;
  70.  
  71.     /* Copy down the sign of y and x */
  72.     signy = 1.;
  73.     if( y < 0 ){y = -y;signy= -1.;}
  74.     if(signx = x < 0 ) x = -x;
  75.  
  76.     if ((t=y/x) < 2.4375) {
  77.         /* Truncate 4(t+1/16) to integer for branching */
  78.         switch (k = t*4 + 0.25) {
  79.  
  80.         case 0:
  81.         case 1: /* T is in [0,7/16] */
  82.             hi = 0;
  83.             lo = 0;
  84.             break;
  85.  
  86.         case 2: /* T is in [7/16,11/16] */
  87.             t = (y * 2 - x) / (x * 2 + y);
  88.             hi = athfhi;
  89.             lo = athflo;
  90.             break;
  91.  
  92.         case 3:
  93.         case 4: /* T is in [11/16,19/16] */
  94.             t = (y - x) / (x + y);
  95.             hi = 51471 / 65536.;
  96.             lo = .00001303156151080961566;
  97.             break;
  98.  
  99.         default: /* T is in [19/16,39/16] */
  100.             t = (y - x * 1.5) / (x + y * 1.5);
  101.             hi = at1fhi;
  102.             lo = at1flo;
  103.             break;
  104.         }
  105.     } /* End of if (t < 2.4375) */
  106.  
  107.     else { /* T >= 2.4375 */
  108. /* 0/0 produces NaN (agrees with IEEE, not UNIX tradition) */
  109.         t = -x / y;
  110.         hi = 51471 / 32768.;
  111.         lo = .00001303156151080961566 * 2;
  112.     }
  113.  
  114.     z = t * t;
  115.     z2 = z * z;
  116.     z = ((lo - t * z *
  117.             (a1 + z * (a2 + z * (a3 + z * (a4 + z * a5))
  118.                     + z2 * z2 * (a6 + z * (a7 + z * (a8
  119.                         + z * (a9 + z * (a10 + z * a11)))))))
  120.         ) + t) + hi;
  121.  
  122.     return ((signx ? PI - z : z) * signy);
  123. }
  124.